function out = compute_elasticities(out, dK_B, setup)
% computes elasticities wrt. a change in K_B
% when computing the elasticities it needs to
% be taken into account that the change in K_B is so distributed
% over the components K_B_M, K_B_R, K_B_C that we have
% 1. Y_B_M = Y_B_R = Y_B_C
% 2. K_B = K_B_M + K_B_R + K_B_C
% This is taken care of by solving the equations specified by the
% sub-function.

    inp_0 = out; % out is a struct, passing using it as in and
                 % output makes sure that those elements which are
                 % aleady part of out are kept and new elements are
                 % added.
    
    gpops_fun = setup.auxdata.gpops_fun;

    x0(1) = inp_0.K_B_M;
    x0(2) = inp_0.K_B_R;
    x0(3) = inp_0.K_B_C;
    
    fun = @(x) equations(x, inp_0, dK_B, gpops_fun);
    
    fprintf('computing elasticities\n')
    sol = fsolve(fun, x0)
    
    inp_1 = inp_0;
    inp_1.K_B_M = sol(1);
    inp_1.K_B_R = sol(2);
    inp_1.K_B_C = sol(3);
    
    out.K_B_0 = inp_0.K_B_M + inp_0.K_B_R + inp_0.K_B_C;
    
    out.dK_B_M = inp_1.K_B_M - inp_0.K_B_M;
    out.dK_B_R = inp_1.K_B_R - inp_0.K_B_R;
    out.dK_B_C = inp_1.K_B_C - inp_0.K_B_C;
    out.dK_B = out.dK_B_M + out.dK_B_R + out.dK_B_C;

    out.Y_M_0 = gpops_fun.prod_fun.Y_M(inp_0);
    out.Y_R_0 = gpops_fun.prod_fun.Y_R(inp_0);
    out.Y_C_0 = gpops_fun.prod_fun.Y_C(inp_0);
    out.Y_B_M_0 = gpops_fun.prod_fun.Y_B_M(inp_0);
    out.Y_B_R_0 = gpops_fun.prod_fun.Y_B_R(inp_0);
    out.Y_B_C_0 = gpops_fun.prod_fun.Y_B_C(inp_0);
    out.Y_E_0 = gpops_fun.prod_fun.Y_E(inp_0);
    out.Y_S_0 = gpops_fun.prod_fun.Y_S(inp_0);
    
    out.Y_M_1 = gpops_fun.prod_fun.Y_M(inp_1);
    out.Y_R_1 = gpops_fun.prod_fun.Y_R(inp_1);
    out.Y_C_1 = gpops_fun.prod_fun.Y_C(inp_1);
    out.Y_B_M_1 = gpops_fun.prod_fun.Y_B_M(inp_1);
    out.Y_B_R_1 = gpops_fun.prod_fun.Y_B_R(inp_1);
    out.Y_B_C_1 = gpops_fun.prod_fun.Y_B_C(inp_1);
    out.Y_E_1 = gpops_fun.prod_fun.Y_E(inp_1);
    out.Y_S_1 = gpops_fun.prod_fun.Y_S(inp_1);    
    
    out.dY_M = out.Y_M_1 - out.Y_M_0;
    out.dY_R = out.Y_R_1 - out.Y_R_0;
    out.dY_C = out.Y_C_1 - out.Y_C_0;
    out.dY_B_M = out.Y_B_M_1 - out.Y_B_M_0;
    out.dY_B_R = out.Y_B_R_1 - out.Y_B_R_0;
    out.dY_B_C = out.Y_B_C_1 - out.Y_B_C_0;
    out.dY_E = out.Y_E_1 - out.Y_E_0;
    out.dY_S = out.Y_S_1 - out.Y_S_0;
    
    out.dY_M_dK_B = out.dY_M/out.dK_B;
    out.dY_R_dK_B = out.dY_R/out.dK_B;
    out.dY_C_dK_B = out.dY_C/out.dK_B;
    out.dY_B_dK_B = out.dY_B_M/out.dK_B;
    out.dY_E_dK_B = out.dY_E/out.dK_B;
    out.dY_S_dK_B = out.dY_S/out.dK_B;
    
    out.eps_Y_M_K_B = out.dY_M/out.dK_B * out.K_B_0/out.Y_M_0;
    out.eps_Y_R_K_B = out.dY_R/out.dK_B * out.K_B_0/out.Y_R_0;
    out.eps_Y_C_K_B = out.dY_C/out.dK_B * out.K_B_0/out.Y_C_0;
    
    out.eps_Y_M_Y_R_K_B = out.eps_Y_M_K_B - out.eps_Y_R_K_B;
    out.eps_Y_C_Y_R_K_B = out.eps_Y_C_K_B - out.eps_Y_R_K_B;
end

function out = equations(x, inp, dK_B, gpops_fun)
    dK_B_M = x(1) - inp.K_B_M;
    dK_B_R = x(2) - inp.K_B_R;
    dK_B_C = x(3) - inp.K_B_C;
    
    inp.K_B_M = x(1);
    inp.K_B_R = x(2);
    inp.K_B_C = x(3);
    
    Y_B_M = gpops_fun.prod_fun.Y_B_M(inp);
    Y_B_R = gpops_fun.prod_fun.Y_B_R(inp);
    Y_B_C = gpops_fun.prod_fun.Y_B_C(inp);
    
    out(1) = Y_B_M - Y_B_R;
    out(2) = Y_B_M - Y_B_C;
    out(3) = dK_B_M + dK_B_R + dK_B_C - dK_B;
end